OPEN Q ACCESS Freely available online 



i»PLOS 



ONE 



Nonlinear Dual Reconstruction of SPECT Activity and 
Attenuation Images 

Huafeng Liu 1 , Min Guo 1 , Zhenghui Hu 1 , Pengcheng Shi 2 , Hongjie Hu 3 * 

1 State Key Laboratory of Modern Optical Instrumentation, Department of Optical Engineering, Zhejiang University, Hangzhou, China, 2B. Thomas Golisano College of 
Computing and Information Sciences, Rochester Institute of Technology, Rochester, New York, United States of America, 3 Department of Radiology, Sir Run Run Shaw 
Hospital, College of Medicine, Zhejiang University, Hangzhou, China 



CrossMark 



Abstract 

In single photon emission computed tomography (SPECT), accurate attenuation maps are needed to perform essential 
attenuation compensation for high quality radioactivity estimation. Formulating the SPECT activity and attenuation 
reconstruction tasks as coupled signal estimation and system parameter identification problems, where the activity 
distribution and the attenuation parameter are treated as random variables with known prior statistics, we present a 
nonlinear dual reconstruction scheme based on the unscented Kalman filtering (UKF) principles. In this effort, the dynamic 
changes of the organ radioactivity distribution are described through state space evolution equations, while the photon- 
counting SPECT projection data are measured through the observation equations. Activity distribution is then estimated 
with sub-optimal fixed attenuation parameters, followed by attenuation map reconstruction given these activity estimates. 
Such coupled estimation processes are iteratively repeated as necessary until convergence. The results obtained from 
Monte Carlo simulated data, physical phantom, and real SPECT scans demonstrate the improved performance of the 
proposed method both from visual inspection of the images and a quantitative evaluation, compared to the widely used 
EM-ML algorithms. The dual estimation framework has the potential to be useful for estimating the attenuation map from 
emission data only and thus benefit the radioactivity reconstruction. 
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Introduction 

Single photon emission computed tomography (SPECT) has 
become an indispensable tool in clinical trials and medical 
practice. Attenuation correction in SPECT has significant values 
to better understand the physiological processes associated with 
the disease (i.e. cancers, heart diseases) and to provide improve- 
ment in patient diagnosis and treatment. However, even with 
ample efforts devoted to the development of various attenuation 
estimation and compensation techniques, it remains one of the key 
open issues in SPECT imaging [1-1 1]. 

Tissue attenuation map is usually estimated based on transmis- 
sion data by scanning the patient with a rotating external 
radionuclide source [3,5,6,12], or obtained from X-ray computed 
tomography (CT) system [10,11,13-16]. However, transmission 
based attenuation correction clearly increases the patient's dose, 
and requires maintaining additional radioactive sources. Further, 
if multiple imaging sessions are needed, it may be difficult for some 
patients to tolerate for a longer scan at one time, and leads to co- 
registration problem in emission image reconstruction, especially 
for deformed tissues and organs. In addition to the added 



equipment cost and the well-known beam hardening problem, 
similar registration issue exists for CT attenuation data. 

It has been the goal of many recent research efforts to 
simultaneously estimate the activity and attenuation distributions 
from emission data. Iterative statistical methods have been 
extensively studied, with the main incentive being that they 
explicitly take into account the specific SPECT data statistics. 
Some of the most notable works include the use of differential 
attenuation method [17], gradient ascent [18], Tikhonov regular- 
ization [19], and expectation maximization (EM) [20—23]. With 
the SPECT imaging and measurement processes in state space 
representation, a recent work has adopted the extended Kalman 
filtering (EKF) procedures to linearize the augmented state 
representation to provide the joint estimates in the minimum- 
mean-square-error sense [24]. Such EKF based framework, 
however, has several potential drawbacks. First, the derivation of 
the Jacobian matrices, the linear approximations to the nonlinear 
functions, often leads to filter instability. Furthermore, it has been 
shown that estimation bias may originate from a coupling between 
the state variables and the model parameters, which suggests that 
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if the system parameters can be separated from the system state 
variables, the precision of the estimation could be improved. 

Following the same spirit while addressing the limitations, we 
present a robust unscented Kalman filter framework for the joint 
estimation of SPECT activity and attenuation parameters. Instead 
of linearizing using the Jacobian matrices and thus overcoming the 
associated shortcomings, our effort deals with the nonlinear 
estimation process by using a deterministic sampling approach 
to capture the mean and covariance estimates [25]. In addition, 
two iteratively coupled filters are used to sequentially estimate the 
activity variables (with fixed attenuation estimates from the last 
iteration) and the attenuation map (with fixed activity values from 
the previous estimation). Furthermore, such framework can 
explicitly recognize the uncertainties in the measurement data 
and the model structure, thus has the potential to produce more 
robust reconstructions. 

Materials and Methods 

SPECT Emission Scan Model 

In SPECT imaging, once radio-tracers are injected into the 
subjects, they are delivered to the tissues/organs by the blood flow 
and participate in the related physiologic/metabolic processes. 
The general mathematical model for the emission scan can be 
stated as 



y(t) = Gx(t) + e 0 



(1) 



where column vector x(t) = {xj\j '= l,...,N} is the radiopharma- 
ceutical concentration of the object, e 0 is the background noises 
(e.g. scatter events), and y is the emission measurement data that is 
acquired by a rotating detector head around the patient at each 
angle and represented lexicographically as a column vector 
y{t) = {yi\i=\,...,M}. Here, i indicates different projection 
defined by rotating angle and different detector bin, and M is 
the total number of projections. G represents the emission system 
matrix including attenuation effects, and its element indicates the 
probability that a photon emitted from pixel gets detected in a 
specific detector bin. 

Let the attenuation map be given by column vector 
fi = {/i k \k= l,...,iV}, we may explicitly account for the attenuation 
effects in the system matrix G by factoring it as 



G = {A'D)=[e-^'D s j 



(2) 



Here, The symbol '•' means the dot product of two matrices. 
\L\ i j k = lijk represents the length of the ray through voxel k with 
respect to a photon emitted from voxel i being detected in 
projection bin j. The photon survival probability considering 
attenuation can now be represented by a matrix A with elements 

\ A \ij = a y as 



(3) 



The term D, with elements [D] i j = dij, is the photon-detection 
probability without taking into account the attenuation effects. 

The mathematical model of the SPECT emission measurement 
(Eqn. (1)) can now be rewritten as 



y=Gx + e 0 =(e~ [Lli] -D > )x+e 0 
or in discrete form 



ex P -^'iikPk ) d v x j 



(4) 



(5) 



In the following sections, we will present the dual estimation 
method which works by alternating between using the UKF state- 
filter to estimate the radioactivity based on fixed attenuation 
parameters, and the UKF parameter-filter to estimate the 
attenuation coefficients given previously estimated activity states. 
Although two separate state-space representations are constructed 
for the activity and attenuation estimation problems, both of them 
use the same emission scan model (Eqn. (4) or (5)) as the 
measurement equation. 

Activity Estimation: UKF State Filter 

State Space Representation for Radioactivity. In emis- 
sion tomography, the goal is to reconstruct the radioactivity 
distribution x(t) from the measurement data y(t). The system 
equation of the SPECT imaging system, which describes the 
radioactivity evolution of the pixels, can be written in the form of 



x(t+l) = Sx(t) + D s 



(6) 



with initial activity Xq and system noise v s that accounts for the 
statistical uncertainty of the imaging model. In general, Eqn. (6) 
represents the dynamic changes of the state variable x(f), and it 
reduces to the conventional static reconstruction problem when 
the transition matrix S is an identity. The associated measurement 
equation, which describes the observations provided by the 
imaging data y(t), is given by Eqn. (4). And now we have the 
following state space representation for radioactivity distribution: 



X(t+I) = x(t) + D s 



y(t) = G(fi(t))x(t) + e 0 = h(x(t),fi(t)) + e 0 



(7) 



(8) 



where v s and e 0 , with covariance matrices Q s and R 0 , model the 
uncertainties of the imaging system and the measurement data 
respectively. 

Two important observations on the nonlinear SPECT imaging 
system represented by Eqns. (7) and (8). First, since the noises in 
emission sinogram are typically Poisson distributed, it is difficult to 
perform standard estimation on such non-Gaussian system. By 
applying the Anscombe transformation [26], however, the Poisson 
noise could be converted into a Gaussian one, and various Kalman 
filtering techniques become viable options. Secondly, the EKF, 
probably the most widely used estimation algorithm for nonlinear 
systems, has several drawbacks. EKF requires the system is almost 
linear on the time scale of the updates, and is difficult to 
implement and to set proper parameters due to the cross-talk 
between state and parameter. To overcome such limitations, we 
have adopted the unscented transformation (UT) to accurately 
propagate mean and covariance information through nonlinear 
transformations [25], with little additional computational cost. 
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UKF State-Filter. For the SPECT imaging system given by 
Eqns. (7) and (8), the purpose of the UKF state-filtering is to search 
for the optimal estimates of radio-tracer variable x(t), given fixed 
attenuation values fi(t). Since detailed discussion of the unscented 
Kalman filtering is certainly beyond the scope of this paper [25], 
we present here a more algorithmic description while ignoring 
some theoretical considerations. 

The unscented transformation is a method for calculating the 
statistics of a random variable. Given a random variable x 
(dimension L) with mean x and covariance P x through a nonlinear 
function y = f(x). To compute the statistics of y, we form a sigma 
point matrix X and their weights according to the following: 

Ao = x 



1. Calculate the sigma point weights as in (9); 

2. Project the state variable x(t) ahead: 

x(t\t-\) = x(t-\) 

3. Project the error covariance P x (t) ahead: 

Px(t\t- 1) = Px(t-l) + Qs 

4. Calculate the sigma points as defined in Eqn. (9 

X(t\t-\) = [x(t\t- \\x(t\t- 1) + 



(10) 



(11) 



\/(L+ X)P x {t\t- 1) ,x(t\t - 1)- 



A' / = x-(v / (^ + A)Px)._ l ,/ = ^+1,...,2L 



y/(L + X)PxW-l)] 



(12) 



it o 



(m)_ 



(9) 



w 0 



(c). 



L + \ 



-(l-a 2 + £),/ = 0 



1 



2{L + \Y 



1,...,2L 



where \=ct 1 (L-\-K) — L. The parameter ae[0,l] determines the 
spread of the sigma points, > 0 is used to incorporate any prior 
knowledge about the distribution of x, and K is a scaling parameter 
and is often set to be zero. Once having the above definitions, the 
UKF state-filter is initialized with x(0) = Xo and covariance matrix 
P X (Q), the state estimates and their error covariance matrices are 
computed sequentially until convergence: 



5. Filtering of the measurement equations: 

y(t\t-\)=h(X(t\t-\),fi(t-\)) 



i = 0 



2L 



-m-mya,t\t-i) 



(13) 



(14) 



(15) 




Figure 1. Simulated and physical phantom used in the experiments. From left to right: activity {left) and attenuation {middle) distributions of 
simulated Zubal phantom, physical imaging phantom with three different material rods inside, {right) 
doi:1 0.1 371/journal.pone.01 06951. gOOl 
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Figure 2. Synthetic Data. Top: Attenuation maps recovered by the EKF {left) and dual UKF {right) frameworks for low count measurements. Bottom: 
Horizontal profiles along the 13 th {left) and 19 th {right) rows of the recovered attenuation maps. 
doi:1 0.1 371 /journal. pone.01 06951 .g002 
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Figure 3. Synthetic Data. From left to right: activity maps recovered by FBP, Elvl-ML, and UKF methods for high count measurements. 
doi:1 0.1 371 /journal. pone.01 06951 .g003 
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Figure 4. Synthetic Data. Top: Horizontal profiles along the 13 th {left), W h {middle), and 23 rd {right) rows of the recovered activity maps for low 
count measurements. Bottom: Horizontal profiles along the 13 th {left), 18 th {middle), and 23 rd {right) rows of the recovered activity maps for high count 
measurements. 

doi:1 0.1 371 /journal.pone.01 06951 .g004 



6. Compute the Kalman Gain: 

2L 
i=0 

-x(t\t-my(i,t\t-\) 
-KAt-\)) T 

K x (t) = P x (t)y(t)P~x(f) 

7. Update the estimate with the measurement: 

x(t) = x(t\ t-\) + K x (t)(y(t)-m- 1 )) 

8. Update the error covariance: 

Px(t) = Px(t\t- 1) - K x (t)P yx{t) K x (t) T 



(16) 



(17) 



(18) 



(19) 



Here, the previously stored information in the prediction step is 
combined with the new information coming from the next 
measurement y(t) and the Kalman gain matrix to refine x(t) 
and P x ( t ) in the correction step. The covariance of the 
measurement error R 0 and system error Q s is assumed to be 
known and set to time-invariant. 



Attenuation Estimation: UKF Parameter Filter 

Once the UKF state-filter converges, it is followed by the 
estimation of a coupled UKF parameter-filter aiming to recover 
the attenuation map of the object being imaged, given the 
estimated radioactivity map. The system equations for the 
parameter-filter are 



y(t) = h(x(t),fi(t)) + e 0 



(20) 
(21) 



Here, v p is the process noise with covariance matrix Q p , and we 
assume that the attenuation parameter vector [i is temporally 



Table 1. RMSE values of estimated activity maps for the synthetic data. 







FBP 


EM 


EKF 


UKF 


low count measurement 


0.6724 


0.5469 


0.4152 


0.3660 


high count measurement 


0.5252 


0.3933 


0.3710 


0.3575 



doi:1 0.1 371/journal.pone.01 06951 .t001 
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Figure 5. Reconstructed activity map of the physical phantom by FBP (top lefts, EM-ML (top right), EKF (bottom left) and \JKF(bottom 
right) (arrows indicate cold areas), and the associated color scale. 

doi:1 0.1 371 /journal. pone.01 06951 .g005 



constant. Initializing the unscented parameter filter with fi(0) = jl 0 
and covariance matrix P^(0), the parameter-filter follows similar 
recursion steps as the state-filter (Eqns. (9)— (17)) until convergence, 
given the sigma point calculation scheme of Eqn. 9. 

The coupled radioactivity state and attenuation parameter 
estimation processes are iteratively repeated as necessary, until 
stable results are achieved. The final optimal estimates then 
become the reconstructed SPECT activity and attenuation maps. 

Results 

Validation with Synthetic Data 

Synthetic Zubal phantom is used to quantitatively evaluate the 
accuracy and robustness of the framework, where a simplified 



human thorax with two lungs is represented by 32 x 32 
attenuation and activity (Fig. 1). Specifically, while both lungs 
have average attenuation coefficient of 0.04/ cm, attenuation of 
the left lung is nonuniform but that of the right lung is uniform. 
SPECT projection data have been generated for a parallel beam 
geometry and 90 views uniformly spaced over 360 degrees. To 
generate realistic data, simulations in our study are performed 
using toolbox GATE [27], which can provide a relatively accurate 
reference for the assessment of new image reconstruction 
algorithms. And we have performed two studies, one with mean 
activities of 50 counts/pixel (low count), and the other with 200 
counts /pixel (high count). 

For these two sets of synthetic projection data, the radioactivity 
maps are reconstructed using four reconstruction methods: FBP 
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Figure 6. Reconstructed attenuation map of the physical phantom by EKF (/eft) and UKF (right), and the associate color scale. 

doi:1 0.1 371 /journal.pone.01 06951 .g006 



[28], EM-ML [29], EKF [24], and dual UKF (The attenuation 
maps are also reconstructed for the EKF and UKF methods). For 
the UKF framework, the state (activity map) is estimated by the 
UKF state-filter, where, for the first iteration, the attenuation 
parameters come from the initial guess. Otherwise, we use the 
values estimated in the last iteration from the UKF parameter- 
filter. Consequently, these state estimates are used by the UKF 
parameter-filter to recover the attenuation coefficients. This 
process is executed iteratively until it meets the convergence 
criterion, which is defined using two consecutive normalized errors 
%(t+ 1) and x(f) through 1) — #(011 <Q with q being a small 

constant, and % defines the normalized error between the 
estimated and the exact value (S = x for the state-filter process, 
H = fi for the parameter-filter process) with 



i Ef=iis;-sn 2 



N 



(22) 



where H • is the estimated value, is the corresponding true value, 
and i indicates the pixel. 

A detailed statistical analysis on the estimation results against 
the ground truth phantom map is performed. Let N p and x be the 
total number of pixels in the region of interest(ROI), e.g. body 
anatomy, and the final reconstruction results respectively, and x tr 
be the ground truth, we have the following error definitions: 



RMSE=[ — J2(x-x lr Y 



0.5 



(23) 



The results for the reconstruction of the attenuation maps are 
shown in Fig. 2 for the low count case. Visually, the UKF- 
reconstructed attenuation image clearly shows all relevant 
anatomical structures, where the lungs are easily seen with clear 
shape and size. The RMSE values in ROI of body anatomy using 
the UKF method are 0.0190 (low counts) and 0.0187 (high 





Figure 7. Reconstructed attenuation (left, by UKF) and activity maps of the patient by UKF (middle) and FBP (right), and the 
associated color scales. 

doi:1 0.1 371 /journal.pone.01 06951 .g007 
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counts), which, though not as impressive, are still somewhat 
smaller than the corresponding values of 0.0204 (low counts) and 
0.0199 (high counts) obtained by the EKF method (note that the 
true average attenuation coefficient is 0.04). The worse perfor- 
mance of the EKF strategy could be caused by the first-order 
Taylor approximations of state transition that provided insuffi- 
ciently accurate representations. 

The results for the reconstruction of the activity maps are 
presented in Fig. 3, and the quantitative tabulation of the 
reconstruction accuracy is in Table 1. These figures and results 
illustrate that traditional EM-ML and FBP methods, with 
unknown attenuation map settings, produce some noticeable 
errors. The dual UKF estimation framework, on the other hand, 
consistently yields the best quality radioactivity estimates for both 
high and low count data. Same conclusion can be drawn from the 
visual examples of the selected horizontal profiles, as shown in 
Fig. 4. 

Reconstruction from Physical Phantom Scanning Data 

The second data set used for validation has been on a real 
cylinder phantom. The dimension of the phantom is 200 mm 
(diameter) x 290 mm (depth). A Teflon rod and two hollow 
PMMA cylinders with diameters of 50mm are inserted in the 
phantom's volume, as shown in Fig. 1 . The phantom is filled with 
99m Tc concentration with a total radioactivity of 20mCi (100kBq/ 
cc) and the two hollow cylinder rods are filled with air and pure 
water respectively. The phantom was scanned with a Siemens 
JLCa.m duet ECT scanner by two detector head rotating at total 64 
angle position around 180 degree and the acquiring time at each 
position is 30 seconds. The final sinogram data has 64x128 
projections for each slice. Once again, FBP, EM-ML, EKF and 
dual UKF strategies have been used to reconstruct activity maps 
from the measurement data, as shown in Fig. 5. Visually, it is 
evident that the UKF method produces the best reconstruction 
results, especially for the three cold areas. The results of EKF and 
UKF reconstructed attenuation map are shown in Fig. 6, with 
RMSE values of 0.0007 for UKF and 0.0044 for EKF. 

Reconstruction from Real Patient Scanning Data 

The dual reconstruction strategy has also been evaluated on 
clinical studies, where the patients are undergoing 99m Tc sestamibi 
stress tests. Using a Siemens ECam^ scanner, all projections are 
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